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Abstract. An exact algorithm is presented for solving edge weighted graph partitioning prob- 
lems. The algorithm is based on a branch and bound method applied to a continuous quadratic 
programming formulation of the problem. Lower bounds are obtained by decomposing the objective 
function into convex and concave parts and replacing the concave part by an affine underestimate. 
It is shown that the best affine underestimate can be expressed in terms of the center and the radius 
of the smallest sphere containing the feasible set. The concave term is obtained cither by a constant 
diagonal shift associated with the smallest eigenvalue of the objective function Hessian, or by a di- 
agonal shift obtained by solving a semidefinite programming problem. Numerical results show that 
the proposed algorithm is competitive with state-of-the-art graph partitioning codes. 
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1. Introduction. Given a graph with edge weights, the graph partitioning prob- 
lem is to partition the vertices into two sets satisfying specified size constraints, while 
minimizing the sum of the weights of the edges that connect the vertices in the two 
sets. Graph partitioning problems arise in many areas including VLSI design, data 
mining, parallel computing, and sparse matrix factorizations [21 [22J [27J • The 
graph partitioning problem is NP-hard [11] . 

There are two general classes of methods for the graph partitioning problem, 
exact methods which compute the optimal partition, and heuristic methods which 
try to quickly compute an approximate solution. Heuristic methods include spectral 
methods [TJ], geometric methods [T2J, multilevel schemes [T5], optimization-based 
methods [5J, and methods that employ randomization techniques such as genetic 
algorithms [34]. Software which implements heuristic methods includes Metis ( [24l 
123 Hi]), Chaco [HJ, Party [32J, PaToH g], SCOTCH [31], Jostle [37], Zoltan [5J, and 
HUND [13]. 

This paper develops an exact algorithm for the graph partitioning problem. In 
earlier work, Brunetta, Conforti, and Rinaldi [3J propose a branch-and-cut scheme 
based on a linear programming relaxation and subsequent cuts based on separation 
techniques. A column generation approach is developed by Johnson, Mehrotra, and 
Nemhauser [25], while Mitchell [25] develops a polyhedral approach. Karisch, Rendl, 
and Clausen [23j develop a branch-and-bound method utilizing a semidefinite pro- 
gramming relaxation to obtain a lower bound. Sensen [33] develops a branch-and- 
bound method based on a lower bound obtained by solving a multicommodity flow 
problem. 
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In this paper, we develop a branch-and-bound algorithm based on a quadratic 
programming (QP) formulation of the graph partitioning problem. The objective 
function of the QP is expressed as the sum of a convex and a concave function. 
We consider two different techniques for making this decomposition, one based on 
eigenvalues and the other based on semidefinite programming. In each case, we give 
an affine underestimate for the concave function, which leads to a tractable lower 
bound in the branch and bound algorithm. 

The paper is organized as follows. In Scction[2]we review the continuous quadratic 
programming formulation of the graph partitioning problem developed in [14] and we 
explain how to associate a solution of the continuous problem with the solution to the 
discrete problem. In Section [3] we discuss approaches for decomposing the objective 
function for the QP into the sum of convex and a concave functions, and in each case, 
we show how to generate an affine lower bound for the concave part. Section 3] gives 
the branch-and-bound algorithm, while Section [5] provides necessary and sufficient 
conditions for a local minimizer. Section [6] compares the performance of the new 
branch-and-bound algorithm to earlier results given in [23) and [33] . 

Notation. Throughout the paper, || ■ || denotes the Euclidian norm. 1 is the 
vector whose entries are all 1. The dimension will be clear from context. If A 6 R nx ™, 
A >z means that A is positive semidefinite. We let denote the i-th column of the 
identity matrix; again, the dimension will be clear from context. If S is a set, then 
|5| is the number of elements in S. The gradient V/(x) is a row vector. 

2. Continuous quadratic programming formulation. Let G be a graph 
with n vertices 



and let ay be a weight associated with the edge (i, j). When there is no edge between 
i and j, we set ay — 0. For each i and j, we assume that an = and ay = dji] in 
other words, we consider an undirected graph without self loops (a simple, undirected 
graph). The sign of the weights is not restricted, and in fact, ay could be negative, as 
it would be in the max-cut problem. Given integers I and u such that < / < u < n, 
we wish to partition the vertices into two disjoint sets, with between I and u vertices 
in one set, while minimizing the sum of the weights associated with edges connecting 
vertices in different sets. The edges connecting the two sets in the partition are referred 
to as the cut edges, and the optimal partition minimizes the sum of the weights of the 
cut edges. Hence, the graph partitioning problem is also called the min-cut problem. 

In [14] we show that for a suitable choice of the diagonal matrix D, the graph 
partitioning problem is equivalent to the following continuous quadratic programming 
problem: 



where A is the matrix with elements ay. Suppose x is binary and let us define the 
sets 



V = {l,2,---,n}, 



(2.1) 



minimize /(x) := (1 — x) T (A + D)x 
subject to < x < 1, I < l T x < u, 



(2.2) 



V = {i : Xi 



0} and Vi = {i : x t = 1}. 



It can be checked that /(x) is the sum of the weights of the cut edges associated 
with the partition (12. 2|) . Hence, if we add the restriction that x is binary, then (|2.ip 



GRAPH PARTITIONING 



3 



is exactly equivalent to finding the partition which minimizes the weight of the cut 
edges. Note, though, that there are no binary constraints in (|2.ip . The equivalence 
between (|2.1[) and the graph partitioning problem is as follows (see [HI Thm. 2.1]): 

Theorem 2.1. If the diagonal matrix D is chosen so that 
(2.3) da + djj > 2<2y and da > 

for each i and j , then ^2.1\ has a binary solution x and the partition given by \2. 2\ 
is a min-cut. 

The generalization of this result to multiset partitioning is given in [15]. The 
condition (|2 .3(1 is satisfied, for example, by the choice 

djj = max {0, ay, 02 j , . . . , a n j} 

for each j . The proof of Theorem 12.11 was based on showing that any solution to (|2.1ll 
could be transformed to a binary solution without changing the objective function 
value. With a modification of this idea, any feasible point can be transformed to a 
binary feasible point without increasing the objective function value. We now give a 
constructive proof of this result, which is used when we solve (|2.ip . 

Corollary 2.2. If x is feasible in 012.1$ and the diagonal matrix D satisfies 
\2. 3\i . then there exists a binary y with /(y) < /(x) and yi = xi whenever Xi is 
binary. 

Proof. We first show how to find z with the property that z is feasible in (12. 1|) . 
/(z) < /(x), l T z is integer, and the only components of z and x which differ are the 
fractional components of x. If l T x = u or l T x = I, then we are done since / and u are 
integers; hence, we assume that I < l T x < u. If all components of x are binary, then 
we are done, so suppose that there exists a nonbinary component Xi. Since an = 0, a 
Taylor expansion of / gives 

/(x + aei) = /(x) + aV/(x), - a 2 du, 

where is the i-th column of the identity matrix. The quadratic term in the expansion 
is nonpositive since da > 0. If the first derivative term is negative, then increase a 
above until either x% +a becomes 1 or l T x+a is an integer. Since the first derivative 
term is negative and a > 0, /(x + aej) < /(x). If l T x + a becomes an integer, then 
we are done. If Xj + a becomes 1, then we reach a point xi with one more binary 
component and with an objective function value no larger than /(x). If the first 
derivative term is nonnegative, then decrease a below until either xi + a becomes 
or l T x + a is an integer. Since the first derivative term is nonnegative and a < 0, 
/(x+aej) < /(x). If l T x+a becomes an integer, then we are done. If Xi + a becomes 
0, then we reach a point xi with one more binary component and with a smaller value 
for the cost function. In this latter case, we choose another nonbinary component of 
Xi and repeat the process. Hence, there is no loss of generality in assuming that l T x 
is an integer. 

Suppose that x is not binary. Since l T x is an integer, x must have at least two 
nonbinary components, say Xi and Xj. Again, expanding / is a Taylor series gives 

/(x + a(e t - ej)) = /(x) + a(V/(x) ( - V/(x) 3 -) + a 2 (2a 4J - d u - d ]3 ). 

By (I2.3|) . the quadratic term is nonpositive for any choice of a. If the first derivative 
term is negative, then we increase a above until either x% +a reaches 1 or Xj — a reach 



4 



W. W. HAGER, D. T. PHAN, H. ZHANG 



0. Since the first derivative term is negative and a > 0, we have /(x + a(e, — ej)) < 
/(x). If the first derivative term is nonnegative, then we decrease a below until 
either Xi + a reaches or Xj — a reach 1. Since the first derivative term is nonnegative 
and a < 0, it follows that /(x + a(ej — ej)) < /(x). In either case, the value of the 
cost function does not increase, and we reach a feasible point xi with l T xi integer 
and with at least one more binary component. If xi is not binary, then xi must have 
at least two nonbinary components; hence, the adjustment process can be continued 
until all the components of x are binary. These adjustments to x do not increase the 
value of the cost function and we only alter the fractional components of x. This 
completes the proof. □ 

3. Convex lower bounds for the objective function. We compute an exact 
solution to the continuous formulation (|2.ip of graph partitioning problem using a 
branch and bound algorithm. The bounding process requires a lower bound for the 
objective function when restricted to the intersection of a box and two half spaces. 
This lower bound is obtained by writing the objective function as the sum of a convex 
and a concave function and by replacing the concave part by the best affine underes- 
timate. Two different strategies are given for decomposing the objective function. 

3.1. Lower bound based on minimum eigenvalue. Let us decompose the 
objective function /(x) = (1 — x) T (A + D)x in the following way: 

/(x) = (/(x)+a||x|| 2 )- ( r||x|| 2 , 

where a is the maximum of and the largest eigenvalue of A + D. This represents 
a DC (difference convex) decomposition (see [20]) since /(x) +<t||x|| 2 and cr||x|| 2 are 
both convex. The concave term — ||x|| 2 is underestimated by an affine function I to 
obtain a convex underestimate Jl of / given by 

(3.1) / L (x) = (/(x) + ( r||x|| 2 )+^(x). 

We now consider the problem of finding the best affine underestimate I for the con- 
cave function — ||x|| 2 over a given compact, convex set denoted C. The set of affine 
underestimators for — ||x|| 2 is given by 

Si = {£ : K" -> R such that I is affine and - ||x|| 2 > £(x) for all x e C}. 

The best affine underestimate is a solution of the problem 

(3.2) mm max - (||x|| 2 + *(x)) . 

ttoi xtt 

The following result generalizes Theorem 3.1 in |I6j where we determine the best 
affine underestimate for — ||x|| 2 over an ellipsoid. 

Theorem 3.1. Let C C M™ be a compact, convex set and let c be the center and r 
be the radius of the smallest sphere containing C. This smallest sphere is unique and 
a solution of l\3.2fy is 

r(x) = -2c T x+||c|| 2 -r 2 . 

Furthermore, 

min max - (||x|| 2 + t (x)) = r 2 . 
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Proof. To begin, we will show that the minimization in (|3.2p can be restricted to a 
compact set. Clearly, when carrying out the minimization in (|3.2[) . we should restrict 
our attention to those £ which touch the function /i(x) := — ||x|| 2 at some point in C. 
Let y G C denote the point of contact. Since h{x) > £(x) and h(y) = £(y), a lower 
bound for the error h(x) — £(x) over x G C is 

h(x) - *(x) > |*(x) - £(y)| - \h(x) - h(y)\. 

If M is the difference between the maximum and minimum value of h over C, then 
we have 

(3.3) fc(x)-*(x)> |£(x)-£(y)|-M. 

An upper bound for the minimum in (|3.2p is obtained by the linear function £q 
which is constant on C, with value equal to the minimum of h(x) over xeC. If w is 
a point where h attains its minimum over C, then we have 

max h(x) — £n(x) = max h(x) — h(w) = M. 
xec xec 

Let us restrict our attention to the linear functions £ which achieve an objective 
function value in (|3.2p which is at least as small as that of £q. For these £ and for 
x G C, we have 

(3.4) h(x) - £(x) < max h(x) - < max h(x) - l (x) = M. 

xec xec 

Combining p.3p and (|3.4p gives 

(3.5) |%)-£(y)| <2M. 



Thus, when we carry out the minimization in (|3 . 2[) . we should restrict our attention 
to linear functions which touch h at some point y G C and with the change in £ across 
C satisfying the bound (|3.5|) for all x G C This tells us that the minimization in (|3.2I) 
can be restricted to a compact set, and that a minimizer must exist. 

Suppose that £ attains the minimum in (|3.2p . Let z be a point in C where 
/i(x) — ^(x) achieves its maximum. A Taylor expansion around x = z gives 

h(x) - £(x) = h{z) - £(z) + (Vh(z) - V£)(x - z) - ||x - z|| 2 . 

Since £ G Si, h(x) - £(x) > for all x G C. It follows that 

(3.6) h(z) ~ £{z) > -(Wh(z) - V£)(x - z) + ||x - zf. 
Since C is convex, the first-order optimality conditions for z give 

(Vft(z) - W)(x-z) < 
for all x G C. It follows from (|3.6[) that 

(3.7) /i(z)-£(z)>||x-z|| 2 

for all x G C. There exists x G C such that ||x — z|| > r or else z would be the center 
of a smaller sphere containing C. Hence, (|3.7|) implies that 



h{z)-£{z) > r 2 



() 
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Fig. 3.1. Suppose c C 



It follows that 

(3.8) max h(x) - £(x) > h(z) - £(z) > r 2 . 

We now observe that for the specific linear function £* given in the statement of 
the theorem, (|3.8(1 becomes an equality, which implies the optimality of £* in (|3.2j) . 
Expand ft. in a Taylor series around x = c to obtain 

ft(x) = -||c|| 2 - 2c T (x - c) - ||x - c|| 2 
= -2c T x+||c|| 2 -||x-c|| 2 . 

Subtract ^*(x) = — 2c T x + ||c|| 2 — r 2 from both sides to obtain 

(3.9) fc(x)-r(x) =r 2 - ||x-c|| 2 . 

If c € C, then the maximum in (13. 9|) over x € C is attained by x = c for which 

h(c)-£*(c) = r 2 . 

Consequently, (|3.8p becomes an equality for £ — £*, which implies the optimality of 
£* in (|X2l) . 

We can show that c G C as follows: Suppose c ^ C. Since C is compact and 
convex, there exists a hyperplane 7i strictly separating c and C - see Figure 13.11 If c' 
is the projection of c onto TL, then 

(3.10) ||x-c'|| < ||x-c|| forallxeC. 

Let x' £ C be the point which is farthest from c' and let x G C be the point farthest 
from c. Hence, ||x — c|| = r. By (|3.10j) . we have ||x' — c'| < ||x — c|| = r; it follows 
that the sphere with center c' and radius ||x' — c'|| contains C and has radius smaller 
than r. This contradicts the assumption that r was the sphere of smallest radius 
containing C. 
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The uniqueness of the smallest sphere containing C is as follows: Suppose that 
there exist two different smallest spheres S\ and S2 containing C. Let £3 be the 
smallest sphere containing Si n 1S2. Since the diameter of the intersection is strictly 
less than the diameter of S\ or 1S2 , we contradict the assumption that S\ and 1S2 were 
spheres of smallest radius containing C. □ 

Remark 1. Although the smallest sphere containing C in Theorem \3.1\ is unique, 
the best linear underestimator o//i(x) = — ||x|| 2 is not unique. For example, suppose 
a and b £ R" and C is the line segment 

C = {xel":x = aa|(l- a)b, ae[0,l]}. 

Along this line segment, h is a concave quadratic in one variable. The best affine 
underestimate along the line segment corresponds to the line connecting the ends of 
the quadratic restricted to the line segment. Hence, in K n+ , any hyperplane which 
contains the points (/i(a),a) and (h(h),h) leads to a best affine underestimate. 
Remark 2. Let C be the box 

B = {x e K" : p < x < q}. 

The diameter of B, the distance between the points in B with greatest separation, is 
||p — q||. Hence, the smallest sphere containing B has radius at least ||p — q||/2. // 
xeS, then 

\xi - (Pi + «i)/2| < (qt -Pi)/2 

for every i. Consequently, ||x — (p + q)/2|| < ||p — q||/2 and the sphere with center 
c = (p + q)/2 and radius r — ||p — q||/2 contains B. It follows that this is the smallest 
sphere containing B since any other sphere must have radius at least ||p — q||/2. 

Remark 3. Finding the smallest sphere containing C may not be easy. However, 
the center and radius of any sphere containing C yields an affine underestimate for 
||x|| 2 over C. That is, if S is a sphere with C C S, then the best affine underestimate 
for — ||x|| 2 over S is also an affine underestimate for — ||x|| 2 over C. 

3.2. Lower bound based on semideflnite programming. A different DC 
decomposition of /(x) = (1 — x) T (A + D)x is the following: 

/(x) = (/(x) + x T Ax) - x T Ax, 

where A is a diagonal matrix with i-th diagonal element A,; > 0. We would like to make 
the second term x T Ax as small as possible while keeping the first term /(x) + x T Ax 
convex. This suggests the following semideflnite programming problem 

minimize Y^™, Xi 

(3.11) 

subject to A - (A + D) y 0, A ^ 0, 

where A is the diagonal of A. If the diagonal of A + D is nonnegative, then the 
inequality A y can be dropped since it is implied by the inequality A — (A + D) y 0. 

As before, we seek the best linear underestimate of the concave function — x T Ax 
over a compact, convex set C. If any of the vanish, then reorder the components 
of x so that x = (y, z) where z corresponds to the components of A^ that vanish. Let 
A + be the principal submatrix of A corresponding to the positive diagonal elements, 
and define the set 



C + = {y : (y, z) e C for some z}. 
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The problem of finding the best linear underestimate for — x T Ax over C is essentially 
equivalent to finding the best linear underestimate for — y T A + y over the C + . Hence, 
there is no loss of generality in assuming that the diagonal of A is strictly positive. 
As a consequence of Theorem 13. 11 we have 

Corollary 3.2. Suppose the diagonal of A is strictly positive and let c be the 
center and r the radius of the unique smallest sphere containing the set 

K l ' 2 C := {A 1 /^ : x e C}. 

The best linear underestimate of — x T Ax over the compact, convex set C is 

f (x) = -2c T A 1 / 2 x+ ||c|| 2 -r 2 . 

Furthermore, 

min max - (x T Ax + f (x)) =r 2 , 

where 

S 2 = {(■: R" -> K such that £ is affine and - x T Ax > £(x) for all x G C}. 

Proof. With the change of variables y = A 1//2 x, an affine function in x is trans- 
formed to an affine function in y and conversely, an affine function in y is transformed 
to an affine function in x. Hence, the problem of finding the best affine underestimate 
for — x T Ax over C is equivalent to the problem of finding the best affine underestimate 
for — ||y|| 2 over A 1 / 2 C. Apply Theorem l3.1l to the transformed problem in y, and then 
transform back to x. □ 

Remark 4. If C is the box {x e E™ : < x < 1},. then A 1/2 C is als o a box to 
which we can apply the observation in Remark^ In particular, we have 

(3.12) c = -A 1 / 2 ! = \\ l/2 and r = ||^V 1 / 2 l||/2 = \\X 1/2 \\/2. 

Hence, ||c|| 2 — r 2 = and we have £*{x) = A T x. 
Remark 5. Let us consider the set 

C = {x e R n : < x < 1, l T x = b}, 

where < b < n. Determining the smallest sphere containing A 1/,2 C may not be 
easy. However, as indicated in Remark any sphere containing A 1 / 2 C yields an 
underestimate for x T Ax. Observe that 

A 1 ' 2 C = {yeff":0<y< A 1 / 2 , y^ 1 / 2 = b}. 

As observed in Remark^ the center c and radius r of the smallest sphere S containing 
the set 

{y £ K" : < y < A 1 / 2 } 

are given in \3.12\) . The intersection of this sphere with the hyperplane y T A -1 / 2 = 
b is a lower dimensional sphere S 1 whose center c' is the projection of c onto the 
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hyperplane. S' contains C since C is contained in both the original sphere S and the 
hyperplane. With a little algebra, we obtain 

By the Pythagorean Theorem, the radius r' of the lower dimensional sphere S' is 



\ 




(b- .5n) 2 

Hence, by Corollary 1 3. SX an underestimate of — x T Ax is given by 

*(*) = -A T x+ (^§t) 1 T x+ IIc'H 2 - (r') 2 . 

Since l T x = b when x G C, it can be shown, after some algebra, that £(x) = — A T x 
(all the constants in the affine function cancel). Hence, the affine underestimate I* 
computed in Remark ^ for the unit box and the affine underestimate I computed in 
this remark for the unit box intersect the hyperplane l T x = b are the same. 

4. Branch and bound algorithm. Since the continuous quadratic program 
(12. 1[) has a binary solution, the branching process in the branch and bound algorithm 
is based on setting variables to or 1 and reducing the problem dimension (we do 
not employ bisections of the feasible region as in [16]). We begin by constructing a 
linear ordering of the vertices of the graph according to an estimate for the difficulty 
in placing the vertex in the partition. For the numerical experiments, the order was 
based on the total weight of the edges connecting a vertex to the adjacent vertices. 
If two vertices v\ and i>2 have weights w\ and W2 respectively, then v\ precedes V2 if 

Wl > W2- 

Let vi, V2, . . ., v n denote the ordered vertices. Level i in the branch and bound 
tree corresponds to setting the Vj-th component of x to the values or 1. Each leaf 
at level i represents a specific selection of and 1 values for the v\ through t»j-th 
components of x. Hence, a leaf at level i has a label of the form 

(4.1) t = (bi, b 2 , h), bj = or 1 for 1 < j < i. 

Corresponding to this leaf, the value of the Vj-th component of x is bj for 1 < j < i. 

Let 7^ denote the branch and bound tree at iteration k and let £(Tk) denote the 
leaves in the tree. Suppose r £ £{Tk) lies at level i in 7^ as in (|4.1[) . Let x T denote the 
vector gotten by removing components Vj, 1 < j < i, from x. The Uj-th component 
of x has the pre- assigned binary value bj for 1 < j ' < i. After taking into account 
these assigned binary values, the quadratic problem reduces to a lower dimensional 
problem in the variable x T of the form 

minimize f T (x T ) 
subject to < x r < 1, l T < l T x r < u T , 

where 

i i 

u T = u — bj and l T = I — bj . 

3=1 3 = 1 



10 



W. W. HAGER, D. T. PHAN, H. ZHANG 



Level 1 (vj) 
Level 2 (v 2 ) 
Level 3 (v 3 ) 




Fig. 4.1. Branch and bound tree 



Using the techniques developed in Section [3l we replace f T by a convex lower bound 
denoted and we consider the convex problem 

minimize (x r ) 

(4.2) 

subject to < x r < 1, l T < l T x r < u T . 

Let M(t) denote the optimal objective function value for (I4.2[) . At iteration k, the 
leaf r £ £(Tk) for which M(r) is smallest is used to branch to the next level. If r has 
the form (|4.ip . then the branching processes generates the two new leaves 

(4.3) (b 1 ,b 2 ,...,b l ,0) and (h, b 2 , . . . , b h 1). 

An illustration involving a 3- level branch and bound tree appears in Figure 14.11 

During the branch and bound process, we must also compute an upper bound 
for the minimal objective function value in (|2.1|) . This upper bound is obtained 
using a heuristic technique based on the gradient projection algorithm and sphere 
approximations to the feasible set. These heuristics for generating an upper bound will 
be described in a separate paper. As pointed out earlier, many heuristic techniques are 
available (for example, Metis (;24, 25, 26 ), Chaco [IT], and Party [32]). An advantage 
of our quadratic programming based heuristic is that we start at the solution to the 
lower bounding problem, a solution which typically has fractional entries and which 
is a feasible starting point for (|2.1| . Consequently, the upper bound is no larger than 
the objective function value associated with the optimal point in the lower-bound 
problem. 

Convex quadratic branch and bound (CQB) 

1. Initialize 7o = and k = 0. Evaluate both a lower bound for the solution to 
(12. ip and an upper denoted Uq. 

2. Choose T fc £ 8 {%) such that M(r k ) = min{A/(r) : r £ £{%)}. If M(r k ) = 
Uk, then stop, an optimal solution of (|2.1j) has been found. 

3. Assuming that r k has the form (14. 1|) . let T k +i be the tree obtained by branch- 
ing at Tfc and adding two new leaves as in (|4.3|) : also see Figure [4~T1 Evaluate 
lower bounds for the quadratic programming problems (|4.2[) associated with 
the two new leaves, and evaluate an improved upper bound, denoted Uk+i, by 
using solutions to the lower bound problems as starting guesses in a descent 
method applied to (|2.1|) . 

4. Replace k by k + 1 and return to step 2. 

Convergence is assured since there are a finite number of binary values for the 
components of x. In the worst case, the branch and bound algorithm will build all 
2 n+1 — 1 nodes of the tree. 
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5. Necessary and sufficient optimality conditions. We use the gradient 
projection algorithm to obtain an upper bound for a solution to (|2.1|) . Since the 
gradient projection algorithm can terminate at a stationary point, we need to be able 
to distinguish between a stationary point and a local minimizer, and at a stationary 
point which is not a local minimizer, we need a fast way to compute a descent direction. 

We begin by stating the first-order optimality conditions. Given a scalar A, define 
the vector 



For any vector fi, Af(fi) is a vector of sets whose «-component is the set jV(/Xj). 
The first-order optimality (Karush-Kuhn- Tucker) conditions associated with a local 
minimizer x of (|2.ip can be written in the following way: For some scalar A, we have 

(5.1) 0<x<l, x e 7V(/x(x, A)), Z<l T x<it, and l T x e M{\). 

The first and third conditions in (|5.ip are the constraints in (|2.ip . while the remain- 
ing two conditions correspond to complementary slackness and stationarity of the 
Lagrangian. 

In [T3] we give a necessary and sufficient optimality conditions for (|2.1[) . which 
we now review. Given any x that is feasible in (|2.ip , let us define the sets 

W(x) = {i: Xl = 1}, £(x) = {i : x t = 0}, and .F(x) = {i : < x 4 < 1}. 

We also introduce subsets Uq and Cq defined by 

W o (x,A)={ieW(x): ft (x,A) = 0} and £ (x, A) = {i G £(x) : m(x, A) = 0}. 

Theorem 5.1. Suppose that I = u and D is chosen so that 

(5.2) da + djj > 2a,ij . 

for all i and j . A necessary and sufficient condition for x to be a local minimizer in 
§2. 11 is that the following all hold: 

(PI) For some X, the first-order conditions \5.1\ are satisfied at x. 

(P2) For each i and j £ ^"(x), we have da + djj = 2a^. 

(P3) Consider the three sets i/o(x, A), £o(x, A), and ^"(x). For each i and j in two 
different sets, we have da + djj = 2ay . 
In treating the situation I < u, an additional condition concerning the dual multi- 
pliers A and fi in the first-order optimality conditions (|5.ip enters into the statement 
of the result: 

(P4) If A = /Ltj(x, 0) = for some i, then da — in any of the following three 
cases: 

(a) I < l T x < u. 

(b) Xi > and l T x = u. 



/x(x, A) = (A + D)l - 2(A + D)x + Al, 



and the set-valued maps TV : R — > 2 R and M : R 



2 




R if v = 
{u} if v > 
{/} if v < 
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(c) Xi < 1 and l T x = I. 
COROLLARY 5.2. Suppose that I < u and D is chosen so that 

(5.3) da + d j3 > Zaij ana da > 

for all i and j . A necessary and sufficient condition for x to be a local minimizer in 
(TOl) is that (P1)-(P4) all hold. 

Based on Theorem 15.11 and Corollary 15.21 we can easily check whether a given 
stationary point is a local minimizer. This is in contrast to the general quadratic 
programming problem for which deciding whether a given point is a local minimizer 
is NP-hard (see [29l [30]). We now observe that when x is a stationary point and 
when any of the conditions (P2)-(P4) are violated, then a descent direction is readily 
available. 

Proposition 5.3. Suppose that x is a stationary point for (TJ0]) and (|5.ff[) holds. 
If either (P2) or (P3) is violated, then d = ej — ej, with an appropriate choice of sign, 
is a descent direction. If I < u, A = = /ii(x, 0), and da > 0, then d = e^, with an 
appropriate choice of sign, is a descent direction in any of the cases (a)-(c) o/(P4). 

Proof. The Lagrangian L associated with (|2.1[) has the form 

(5.4) L(x) = /(x) + A(l T x -b)-J2 IH*i - tofa - 1), 

iec ieu 

where b — u if A > 0, b — I if A < 0, and /i. stands for /i.(x, A). The sets C and U 
denote £(x) and U{x.) respectively. By the first-order optimality conditions (|5.ip . we 
have i(x) = /(x) and VL(x) = 0. Expanding the Lagrangian around x gives 

L(x + y) = L(x) + VL(x)y + iy T V 2 L(x)y - /(x) - y T (A + D)y. 
We substitute for L using (|5.4[) to obtain 

/(x + y) = L(x + y) - A(l T (x + y) - b) + ^ mfc + y 4 ) + ^ ^(^ + |fc - 1) 

(5.5) = /(x) - Al T y - y T (A + D)y + ^ Ml2/l + £ 

If (P2) is violated, then there are indices i and j g ^"(x) such that d„ + djj > 2ay. 
We insert y = a(ej — e^) in (|5.5p to obtain 

(5.6) /(x + a(e l - ej)) = /(x) + a 2 (2a tj - d u - djj). 

Since the coefficient of a 2 is negative, d = — ej is a descent direction for the 
objective function. Since < Xi < 1 and < xj < 1, feasibility is preserved for a 
sufficiently small. In a similar manner, if (P3) is violated by indices i and j, then 
(|5.6|) again holds and d = ±(e, — ej) is again a descent direction where the sign is 
chosen appropriately to preserve feasibility. For example, if i £ £o(x) and j E t/o(x), 
then Xi = and Xj = 1. Consequently, x + a(e.; — ej) is feasible if a > is sufficiently 
small. 
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Finally, suppose that I < u, A = = /^(x, 0), and da > 0. Substituting y = ae, 
in (|5.5p yields 

/(x + aej) = /(x) - a 2 d i4 . 

Since the coefficient da of a 2 is positive, d = ±e^ is a descent direction. Moreover, in 
any of the cases (a)-(c) of (P4), x + ad is feasible for some a > with an appropriate 
choice of the sign of d. □ 

We now give a necessary and sufficient condition for a local minimizer to be strict. 
When a local minimizer is not strict, it may be possible to move to a neighboring point 
which has the same objective function value but which is not a local minimizer. 

Corollary 5.4. If x is a local minimizer for (jli.ip and (|5. 3\ holds, then x is a 
strict local minimizer if and only if the following conditions hold: 
(CI) ^"(x) is empty. 

(C2) V/(x)j > V/(x)j for every i G £(x) and j G U(x). 

(C3) If I < u, the first-order optimality conditions (|5. J|) hold for A = 0, and 
Z :={i: V/(x). i = 0} ^ 0, then either 

(a) l T x = u and Xi = for all i G Z or 

(b) l T x = I and xi — 1 for all i G Z. 

Proof. Throughout the proof, we let /x, C and U denote /x(x, A), ^(x), £(x), 
and W(x) respectively, where x is a local minimizer for (|2.1[) and the pair (x, A) satisfies 
the first-order optimality conditions (|5.ip . To begin, suppose that x is a strict local 
minimizer of (|2.1| . That is, /(y) > /(x) when y is a feasible point near x. If T has 
at least two elements, then by (P2) of Theorem 15.11 da + djj = 2ajj for each i and 
j G T . Since the first-order optimality conditions ()5.ip hold at x, it follows from (|5.6|) 
that 

(5.7) /(x + a(e i -e,))=/(x) 

for all a. Since this violates the assumption that x is a strict local minimizer, we 
conclude that < 1. If l T x = u or l T x = I, then since u and / are integers, it is 
not possible for x to have just one fractional component. Consequently, T is empty. 
If / < l T x < u, then by complementary slackness, A = 0. Suppose that \F\ — 1 and 
i G T. By (P4) of Corollary [El d l% = 0. Again, by (JES]) it follows that 

/(x + ae,) = /(x) 

for all a. This violates the assumption that x is a strict local minimizer of (|2.ip . 
Hence, T is empty. 

By the first-order conditions (|5.ip . there exists A such that 

(5.8) m(it, A) > > /ij(x, A) 

for all i G £ and j <E U. If this inequality becomes an equality for some i £ £ and 
.7 G then /i^ = = /Zj-, and by (P3) of Corollary 15. 2[ we have da + djj = 2djj. 
Again, (|5.7p violates the assumption that x is a strict local minimizer. Hence, one of 
the inequalities in (|5.8p is strict. The A on each side of (|5.8p is cancelled to obtain 
(C2). 

Suppose that I < u, X = 0, and Z := {i : V/(x) 4 = 0} ^ 0. When A = 0, we 
have /i(x,0) = V/(x). Hence, Z = {i : ^(x,0) = 0} ^ 0. It follows from (P4) that 
in any of the cases (a)-(c), we have da = 0. In particular, if I < l T x < u, then by 
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(|5.5[) . we have /(x + aej) = /(x) for all a. Again, this violates the assumption that x 
is a strict local minimum. Similarly, if for some i E Z, either Xi > and l T x — u or 
Xi < 1 and l T x = I, the identity /(x + ae{) = /(x) implies that we violate the strict 
local optimality of x. This establishes (C3). 

Conversely, suppose that x is a local minimizer and (C1)-(C3) hold. We will 
show that 

(5.9) V/(x)y > whenever y ^ and x + y feasible in (|2.1[) . 

As a result, by the mean value theorem, /(x + y) > /(x) when y is sufficiently small. 
Hence, x is a strict local minimizer. 

When x + y is feasible in (|2.ip , we have 

(5.10) y l > for alH G £ and y l < for all i e U. 

By the first-order optimality condition (|5.ip . /ij > for all i S £ and (ii < for all 
i ElA. Hence, we have 

(5.11) (V/(x) + Al T )y = /x T y = ]T MiJ/i + ^ WW ^ °- 
We now consider three cases. 

First, suppose that l T y = and y ^0. By (CI) T is empty and hence, by (|5. 10[) . 
Ui > for some i E C and yj < for some j E IA. After adding A to each side in the 
inequality in (C2), it follows that either 

(5.12) min fit > > max fij 
or 



(5.13) min > > max/A,-. 

Combining (|5.1ip . (|5.12[) . and (|5.13|) gives V/(x)y > /Ujj/i — /i^yj > since either 
Hi > or /Zj < 0, and yi > > yj. 

Second, suppose that l T y ^ and A ^ 0. To be specific, suppose that A > 0. By 
complementary slackness, 1 x = u. Since x + y is feasible in flUJ and l T y ^ 0, we 
must have l T y < 0. Hence, by (|5.11|) , V/(x)y > 0. The case A < is similar. 

Finally, consider the case l T y ^ and A = 0. In this case, we must have I < u. If 
the set Z in (C3) is empty, then V/(x)j = m 7^ f° r au h an d by (|5 . llj) . V/(x)y > 0. 
If Z ^ 0, then by (C3), either l T x = u and x t — for all i E Z or l T x = I and Xj = 1 
for all i £ Z. To be specific, suppose that l T x = u and x^ = for all i E Z. Again, 
since x + y is feasible in (|2.1[) and l T y ^ 0, we have l T y < 0. If U = 0, then x = 
since T = %. Since l T y < 0, we contradict the feasibility of x + y. Hence, W / f). 
Since l T y < 0, there exists j E U such that yj < 0. Since Z C C, it follows from 
(|5~T2]) that /Ltj < 0. By ([5TT|) V/(x)y > fj^Uj > °- The case lTx = ' and x i = 1 
for all z E Z is similar. This completes the proof of (15. 9p . and the corollary has been 
established. □ 

6. Numerical results. We investigate the performance of the branch and bound 
algorithm based on the lower bounds in Section [3] using a series of test problems. The 
codes were written in C and the experiments were conducted on an Intel Xeon Quad- 
Core X5355 2.66 GHz computer using the Linux operating system. Only one of the 4 
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processors was used in the experiments. To evaluate the lower bound, we solve (|4.2p 
by the gradient projection method with an exact linesearch and Barzilai-Borwein 
steplength [Tj. The stopping criterion in our experiments was 

||P(x fc -g fc )-x fc || <10" 4 , 

where P denotes the projection onto the feasible set and gk is the gradient of the 
objective function at x&. The solution of the scmidefinite programming problem 
(|3.11|) was obtained using Version 6.0.1 of the CSDP code available at 

https://projects.coin-or.org/Csdp/ 

We compare the performance of our algorithm with results reported by Karisch, Rendl, 
and Clausen in [23] and by Sensen in [33]. Since these earlier results were obtained 
on different computers, we obtained estimates for the corresponding running time on 
our computer using the LINPACK benchmarks [7]. Since our computer is roughly 30 
times faster than the HP 9000/735 used in [53] and it is roughly 7 times faster than the 
Sun UltrSPARC-II 400Mhz machine used in [33], the earlier CPU times were divided 
by 30 and 7 respectively to obtain the estimated running time on our computer. Note 
that the same interior-point algorithm that we use, which is the main routine in the 
CSDP code, was used to solve the semidefmite relaxation in [23] . 

The test problems were based on the graph bisection problem where I — u = n/2. 
Two different data sets were used for the A matrices in the numerical experiments. 
Most of the test problems came from the library of Brunetta, Conforti, and Rinaldi 
[3] which is available at 

ftp://ftp.math.unipd.it/pub/Misc/equicut. 

Some of the test matrices were from the UF Sparse Matrix Library maintained by 
Timothy Davis: 

http://www.cise.ufl.edu/research/sparse/matrices/ 

Since this second set of matrices is not directly connected with graph partitioning, we 
create an A for graph partitioning as follows: If the matrix S from the library was 
symmetric, then A was the adjacency matrix defined as follows: the diagonal of A is 
zero, <Zjj = 1 if Sjj ^ 0, and = otherwise. If S was not symmetric, then A was 
the adjacency matrix of S T S. 

6.1. Lower bound comparison. Our numerical study begins with a compari- 
son of the lower bound of Section I3TT1 based on the minimum eigenvalue of A + D and 
the best affine underestimate, and the lower bound of Section [3721 based on semidef- 
inite programming. We label these two lower bounds LB± and LBi respectively. 
In Table 16.11 the first 5 graphs correspond to matrices from the UF Sparse Matrix 
Library, while the next 5 graphs were from the test set of Brunetta, Conforti, and 
Rinaldi. The column labeled "Opt" is the minimum cut and while n is the prob- 
lem dimension. The numerical results indicate that the lower bound LB2 based on 
semidefmite programming is generally better (larger) than LB\. In Table |6~T1 the best 
lower bound is highlighted in bold. Based on these results, we use the semidefinite 
programming-based lower bound in the numerical experiments which follow. 

6.2. Algorithm performance. Unless stated otherwise, the remaining test 
problems came from the library of Brunetta, Conforti, and Rinaldi [3J. Table W72\ 
gives results for matrices associated with the finite element method [35] ■ The three 
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Table 6.1 
Comparison of two lower bounds 



Graph 


n 


LBi 


LB 2 


Opt 


TinaJDiscal 


11 


0.31 


0.86 


12 


jgl009 


9 


1.55 


1.72 


16 


jglOll 


11 


1.48 


0.94 


24 


Stranke94 


10 


1.76 


1.77 


24 


Hamrlel 


32 


-1.93 


1.12 


17 


4x5t 


20 


-21.71 


5.43 


28 


8x5t 


40 


-16.16 


2.91 


33 


t050 


30 


0.90 


18.54 


397 


2x1 7m 


34 


1.33 


1.27 


316 


s090 


60 


-9.84 


13.10 


238 



methods are labeled CQB (our convex quadratic branch and bound algorithm) , KRC 
(algorithm of Karisch, Rendl, and Clausen [53]), and SEN (algorithm of Sensen [33]). 
"n" is the problem dimension, "%" is the percent of nonzeros in the matrix, and 

nodes" is the number of nodes in the branch and bound tree. The CPU time is 
given in seconds. The best time is highlighted in bold. As can be seen in Table 16.21 
CQB was fastest in 6 out of the 10 problems even though the number of nodes in the 
branch and bound tree was much larger. Thus both KRC and SEN provided much 
tighter relaxations, however, the time to solve their relaxed problems was much larger 
than the time to optimize our convex quadratics. 

Table l6~3l gives results for compiler design problems [TO] [21]. For this test set, 
KRC was fastest in 3 out of 5 test problems. Note though that the times for CQB 
were competitive with KRC. 

Table 16.41 gives results for binary de Bruijn graphs which arise in applications 
related to parallel computer architecture [5] [9] . These graphs are constructed by the 
following procedure. We first build a directed graph using the Mathematica command: 

A = TableForm[ToAdjacencyMatrix[DeBruijnGraph[2, n] ] ] 

To obtain the graph partitioning test problem, we add the Mathematica generated 
matrix to its transpose and set the diagonal to 0. For this test set, SEN had by far 
the best performance. 

Table 1631 gives results for toroidal grid graphs. These graphs are connected with 
an h x k grid, the number of vertices in the graph is n — hk and there are 2hk edges 
whose weights are chosen from a uniform distribution on the interval [1,10]. Since 
Sensen did not solve either this test set, or the remaining test sets, we now compare 
between CQB and KRC. We see in Table ET51 that CQB was faster than KRC in 9 of 
the 10 toroidal grid cases. 

Table RTBl gives results for mixed grid graphs. These are complete graphs associ- 
ated with an planar h x k planar grid; the edges in the planar grid received integer 
weights uniformly drawn from [1,100], while all the other edges needed to complete the 
graph received integer weights uniformly drawn from [1,10]. For these graphs, KRC 
was much faster than CQB. Notice that the graphs in this test set are completely 
dense. One trend that is seen in these numerical experiments is that as the graph 
density increases, the performance of CQB relative to the other methods degrades. 

Results for planar grid graph are given in Table 16.71 These graphs are associated 
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with an h x k grid. There are hk vertices and 2hk — h — k edges whose weights are 
integers uniformly drawn from [1,10]. For this relatively sparse test set, CQB was 
faster in 7 out of 10 problems. 

Table 16.81 gives results for randomly generated graphs. For these graphs, the 
density is first fixed and then the edges are assigned integer weights uniformly drawn 
from [1,10]. For this test set, CQB is fastest in 11 of 20 cases. Again, observe that 
the relative performance of CQB degrades as the density increases, mainly due to the 
large number of nodes in the branch and bound tree. 

7. Conclusions. An exact algorithm is presented for solving the graph parti- 
tioning problem with upper and lower bounds on the size of each set in the partition. 
The algorithm is based on a continuous quadratic programming formulation of the 
discrete partitioning problem. Wc show how to transform a feasible x for the graph 
partitioning QP (|2.1[) to a binary feasible point y with an objective function value 
which satisfies /(y) < /(x). The binary feasible point corresponds to a partition of 
the graph vertices and /(y) is the weight of the cut edges. At any stationary point 
of (|2.ip which is not a local minimizer, Proposition 15.31 provides a descent direction 
that can be used to strictly improve the objective function value. 

In the branch and bound algorithm, the objective function is decomposed into 
the sum of a convex and a concave part. A lower bound for the objective function 
is achieved by replacing the concave part by an affine underestimate. Two different 
decompositions were considered, one based on the minimum eigenvalue of the matrix 
in the objective function, and the other based on the solution to a semidefinite pro- 
gramming problem. The semidefinite programming approach generally led to much 
tighter lower bounds. In a series of numerical experiments, the new algorithm CQB 
(convex quadratic branch and bound) was competitive with state-of-the-art partition- 
ing methods; the relative performance of CQB was better for sparse graphs than for 
dense graphs. 
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Table 6.2 
Mesh Instances 









CQB 




KRC 


SEN 




graph 


n 


% 


#nodes 


time 


#nodes 


time 


#nodes 


time 


m4 


32 


10 


22 


0.05 


1 


0.03 


1 


0.14 


ma 


54 


5 


8 


0.16 


1 


0.10 


1 


0.28 


me 


60 


5 


13 


0.20 


1 


0.13 


1 


0.28 


m6 


70 


5 


205 


0.47 


1 


1.23 


1 


1.43 


mb 


74 


4 


95 


0.43 


1 


0.98 


1 


1.14 


mc 


74 


5 


412 


0.52 


1 


1.53 


1 


1.43 


rnd 


80 


4 


101 


0.55 


1 


0.96 


1 


1.28 


mf 


90 


4 


99 


0.79 


1 


0.80 


1 


1.85 


ml 


100 


3 


200 


1.04 


15 


36.50 


1 


3.00 


m8 


148 


2 


3516 


6.62 


1 


10.70 


1 


4.14 



Table 6.3 
Compiler Design 









CQB 


KRC 


SEN 


graph 


n 


% 


#nodcs time 


#nodes time 


#nodcs time 


cd30 


30 


13 


11 0.05 


1 0.03 


1 0.00 


cd45 


45 


10 


35 0.27 


1 0.23 


1 0.57 


cd47a 


47 


9 


45 0.34 


1 0.33 


7 1.00 


cd47b 


47 


9 


67 0.29 


35 3.73 


3 1.43 


cd61 


61 


10 


95 0.86 


1 0.67 


6 6.00 



Table 6.4 
de Bruijn Networks 









CQB 


KRC 


SEN 


graph 


n 


% 


#nodes time 


#nodes time 


#nodes time 


dcbr5 


32 


12 


57 0.11 


3 0.20 


1 0.00 


debr6 


64 


6 


7327 2.25 


55 15.63 


1 1.00 


dcbr7 


128 


3 


16140945 1:22:45 


711 46:36 


1 10.28 



Table 6.5 

Toroidal Grid: a weighted hxk grid with hk vertices and 2hk edges that received integer weights 
uniformly drawn from [1,10] 









CQB 


KRC 


graph 


n 


% 


#nodes time 


#nodes time 


4x5t 


20 


21 


13 0.01 


1 0.03 


6x5t 


30 


14 


46 0.05 


1 0.10 


8x5t 


40 


10 


141 0.16 


1 0.20 


21x2t 


42 


10 


18 0.02 


1 0.17 


23x2t 


46 


9 


78 0.15 


33 4.16 


4xl2t 


48 


9 


69 0.17 


3 0.56 


5xl0t 


50 


8 


129 0.24 


1 0.20 


6xl0t 


60 


7 


992 0.54 


43 11.66 


7xl0t 


70 


6 


844 0.68 


47 19.06 


10x8t 


80 


5 


420 0.91 


45 31.46 
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Table 6.6 
Mixed Grid Graphs 









CQB 


KRC 


graph 


n 


% 


#nodes time 


#nodes time 


2x10m 


20 


100 


150 0.03 


1 0.03 


6x5m 


30 


100 


2476 0.20 


1 0.03 


2x1 7m 


34 


100 


42410 2.12 


21 0.96 


10x4m 


40 


100 


51713 3.74 


2 0.06 


5x10m 


50 


100 


3588797 296.19 


1 0.06 



Table 6.7 
Planar Grid 









CQB 


KRC 


graph 


n 


% 


#nodcs time 


#nodes time 


10x2g 


20 


15 


10 0.01 


1 0.03 


5x6g 


30 


11 


44 0.05 


1 0.10 


2xl6g 


32 


9 


23 0.06 


1 0.13 


18x2g 


36 


8 


19 0.08 


1 0.06 


2xl9g 


38 


8 


53 0.29 


49 1.83 


5x8g 


40 


9 


24 0.08 


1 0.06 


3xl4g 


42 


8 


31 0.14 


5 0.60 


5xl0g 


50 


7 


178 0.34 


1 0.30 


6xl0g 


60 


6 


224 0.35 


57 10.63 


7xl0g 


70 


5 


271 0.63 


61 18.56 



Table 6.8 
Randomly Generated Graphs 









CQB 




KRC 




graph 


n 


% 


#nodes 


time 


#nodes 


time 


v090 


20 


10 


12 


0.01 


1 


0.03 


vOOO 


20 


100 


952 


0.02 


1 


0.03 


t090 


30 


10 


10 


0.05 


1 


0.03 


t050 


30 


50 


5081 


0.32 


17 


0.73 


tooo 


30 


100 


122670 


3.79 


3 


0.20 


q090 


40 


10 


89 


0.14 


1 


0.13 


q080 


40 


20 


914 


0.24 


31 


2.30 


q030 


40 


70 


554652 


32.23 


23 


2.06 


q020 


40 


80 


1364517 


72.58 


7 


0.83 


qOlO 


40 


90 


4344123 


217.16 


13 


1.36 


qOOO 


40 


100 


8186984 


380.72 


1 


0.13 


c090 


50 


10 


397 


0.29 


1 


0.33 


c080 


50 


20 


14290 


2.20 


45 


6.13 


c070 


50 


30 


136290 


15.70 


49 


8.06 


c030 


50 


70 


22858729 2756.26 


51 


5.46 


c290 


52 


10 


340 


0.34 


1 


0.40 


c490 


54 


10 


1443 


0.54 


15 


3.30 


c690 


56 


10 


3405 


0.82 


3 


1.00 


c890 


58 


10 


13385 


2.66 


71 


17.53 


s090 


60 


10 


8283 


2.01 


37 


9.90 



